{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Quantitative Finance\n",
    "\n",
    "Copyright (c) 2019 Python Charmers Pty Ltd, Australia, <https://pythoncharmers.com>. All rights reserved.\n",
    "\n",
    "<img src=\"img/python_charmers_logo.png\" width=\"300\" alt=\"Python Charmers Logo\">\n",
    "\n",
    "Published under the Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) license. See `LICENSE.md` for details.\n",
    "\n",
    "Sponsored by Tibra Global Services, <https://tibra.com>\n",
    "\n",
    "<img src=\"img/tibra_logo.png\" width=\"300\" alt=\"Tibra Logo\">\n",
    "\n",
    "\n",
    "## Module 1.3: Ordinary Least Squares\n",
    "\n",
    "### 1.3.2 Regression Tests\n",
    "\n",
    "In this module we will look further into Multivariate OLS and examine some of the requirements of the algorithm, as well as some of the details of the regression results we saw in the last module.\n",
    "\n",
    "When performing OLS for Linear Regression Models, there are a few assumptions that need to be met. The key ones are:\n",
    "\n",
    "The first assumption is the key one - that is that the relationship between $X$ and $Y$ can, in fact, be described using the model $Y = X\\beta + u$. It may *not* be able to be precisely modeled this way, but it may be possible to get close enough that it doesn't matter.\n",
    "\n",
    "The second assumption is that the expected value of $u$ is zero. There may be fluctuations in the vector $u$, but the overall expected value is 0. More formally, we assume that $E(u|X) = 0$, that is the expected value of $u$ when given $X$ is zero. If it were not, then we can alter the bias term to make it zero, which would be learned from the OLS, giving us our zero value!\n",
    "\n",
    "The third assumption is that the error term ($u$) and the data itself $X$ do not have any correlation. In other words, $u$ is unexplained error that cannot be explained by the data. Put more formally, there is no heteroskedasticity or autocorrelation between $u$ and $X$, which is a stronger assumption than the second, but along the same lines. We will cover these terms in a later module more formally.\n",
    "\n",
    "The fourth assumption is that $X$ has a finite variance. This is sometimes (slightly incorrectly) referred to as $X$ being non-stochastic. We will investigate how variance plays into the model in several later modules.\n",
    "\n",
    "The fifth assumption is that there are no linear relation between the measurements (variables, columns, features) in $X$, known as having **full column rank**.\n",
    "\n",
    "If any of these assumptions are untrue, the resulting model does not necessarily have the properties we will discuss in the rest of this module, and the model itself might be biased or inaccurate. However, it may still be *useful* in a practical sense. For instance, if two variables are slightly linearly related, we break the last assumption, however in practice the model is generally still useful. However if they are heavily related, then the resulting model will be unstable.\n",
    "\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "    Like most models and concepts, there is always some debate about the definitions and assumptions behind them. Further, some people use the same term to describe different concepts. When discussing an algorithm, it would be best practice to note any key assumptions or variance from the \"norm\" that you consider. If you aren't sure, provide a reference.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run setup.ipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load in some data for a regression problem and have a look at the results. In this dataset, we are trying to predict house prices from other characteristics of the area, in Boston, Massachusetts. Prices are in thousands, but are from 1978, so are quite low!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "\n`load_boston` has been removed from scikit-learn since version 1.2.\n\nThe Boston housing prices dataset has an ethical problem: as\ninvestigated in [1], the authors of this dataset engineered a\nnon-invertible variable \"B\" assuming that racial self-segregation had a\npositive impact on house prices [2]. Furthermore the goal of the\nresearch that led to the creation of this dataset was to study the\nimpact of air quality but it did not give adequate demonstration of the\nvalidity of this assumption.\n\nThe scikit-learn maintainers therefore strongly discourage the use of\nthis dataset unless the purpose of the code is to study and educate\nabout ethical issues in data science and machine learning.\n\nIn this special case, you can fetch the dataset from the original\nsource::\n\n    import pandas as pd\n    import numpy as np\n\n    data_url = \"http://lib.stat.cmu.edu/datasets/boston\"\n    raw_df = pd.read_csv(data_url, sep=\"\\s+\", skiprows=22, header=None)\n    data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])\n    target = raw_df.values[1::2, 2]\n\nAlternative datasets include the California housing dataset and the\nAmes housing dataset. You can load the datasets as follows::\n\n    from sklearn.datasets import fetch_california_housing\n    housing = fetch_california_housing()\n\nfor the California housing dataset and::\n\n    from sklearn.datasets import fetch_openml\n    housing = fetch_openml(name=\"house_prices\", as_frame=True)\n\nfor the Ames housing dataset.\n\n[1] M Carlisle.\n\"Racist data destruction?\"\n<https://medium.com/@docintangible/racist-data-destruction-113e3eff54a8>\n\n[2] Harrison Jr, David, and Daniel L. Rubinfeld.\n\"Hedonic housing prices and the demand for clean air.\"\nJournal of environmental economics and management 5.1 (1978): 81-102.\n<https://www.researchgate.net/publication/4974606_Hedonic_housing_prices_and_the_demand_for_clean_air>\n",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[4], line 3\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;66;03m# Let's load a dataset from the scikit learn repository\u001b[39;00m\n\u001b[0;32m      2\u001b[0m \u001b[38;5;66;03m# scikit-learn is a machine learning library, and has a few sample datasets \u001b[39;00m\n\u001b[1;32m----> 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01msklearn\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mdatasets\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m load_boston\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\sklearn\\datasets\\__init__.py:157\u001b[0m, in \u001b[0;36m__getattr__\u001b[1;34m(name)\u001b[0m\n\u001b[0;32m    108\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m name \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mload_boston\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m    109\u001b[0m     msg \u001b[38;5;241m=\u001b[39m textwrap\u001b[38;5;241m.\u001b[39mdedent(\u001b[38;5;124m\"\"\"\u001b[39m\n\u001b[0;32m    110\u001b[0m \u001b[38;5;124m        `load_boston` has been removed from scikit-learn since version 1.2.\u001b[39m\n\u001b[0;32m    111\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    155\u001b[0m \u001b[38;5;124m        <https://www.researchgate.net/publication/4974606_Hedonic_housing_prices_and_the_demand_for_clean_air>\u001b[39m\n\u001b[0;32m    156\u001b[0m \u001b[38;5;124m        \u001b[39m\u001b[38;5;124m\"\"\"\u001b[39m)\n\u001b[1;32m--> 157\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mImportError\u001b[39;00m(msg)\n\u001b[0;32m    158\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m    159\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mglobals\u001b[39m()[name]\n",
      "\u001b[1;31mImportError\u001b[0m: \n`load_boston` has been removed from scikit-learn since version 1.2.\n\nThe Boston housing prices dataset has an ethical problem: as\ninvestigated in [1], the authors of this dataset engineered a\nnon-invertible variable \"B\" assuming that racial self-segregation had a\npositive impact on house prices [2]. Furthermore the goal of the\nresearch that led to the creation of this dataset was to study the\nimpact of air quality but it did not give adequate demonstration of the\nvalidity of this assumption.\n\nThe scikit-learn maintainers therefore strongly discourage the use of\nthis dataset unless the purpose of the code is to study and educate\nabout ethical issues in data science and machine learning.\n\nIn this special case, you can fetch the dataset from the original\nsource::\n\n    import pandas as pd\n    import numpy as np\n\n    data_url = \"http://lib.stat.cmu.edu/datasets/boston\"\n    raw_df = pd.read_csv(data_url, sep=\"\\s+\", skiprows=22, header=None)\n    data = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])\n    target = raw_df.values[1::2, 2]\n\nAlternative datasets include the California housing dataset and the\nAmes housing dataset. You can load the datasets as follows::\n\n    from sklearn.datasets import fetch_california_housing\n    housing = fetch_california_housing()\n\nfor the California housing dataset and::\n\n    from sklearn.datasets import fetch_openml\n    housing = fetch_openml(name=\"house_prices\", as_frame=True)\n\nfor the Ames housing dataset.\n\n[1] M Carlisle.\n\"Racist data destruction?\"\n<https://medium.com/@docintangible/racist-data-destruction-113e3eff54a8>\n\n[2] Harrison Jr, David, and Daniel L. Rubinfeld.\n\"Hedonic housing prices and the demand for clean air.\"\nJournal of environmental economics and management 5.1 (1978): 81-102.\n<https://www.researchgate.net/publication/4974606_Hedonic_housing_prices_and_the_demand_for_clean_air>\n"
     ]
    }
   ],
   "source": [
    "# Let's load a dataset from the scikit learn repository\n",
    "# scikit-learn is a machine learning library, and has a few sample datasets \n",
    "from sklearn.datasets import load_boston"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston_data = load_boston()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sklearn.utils.Bunch"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(boston_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".. _boston_dataset:\n",
      "\n",
      "Boston house prices dataset\n",
      "---------------------------\n",
      "\n",
      "**Data Set Characteristics:**  \n",
      "\n",
      "    :Number of Instances: 506 \n",
      "\n",
      "    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.\n",
      "\n",
      "    :Attribute Information (in order):\n",
      "        - CRIM     per capita crime rate by town\n",
      "        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.\n",
      "        - INDUS    proportion of non-retail business acres per town\n",
      "        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)\n",
      "        - NOX      nitric oxides concentration (parts per 10 million)\n",
      "        - RM       average number of rooms per dwelling\n",
      "        - AGE      proportion of owner-occupied units built prior to 1940\n",
      "        - DIS      weighted distances to five Boston employment centres\n",
      "        - RAD      index of accessibility to radial highways\n",
      "        - TAX      full-value property-tax rate per $10,000\n",
      "        - PTRATIO  pupil-teacher ratio by town\n",
      "        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town\n",
      "        - LSTAT    % lower status of the population\n",
      "        - MEDV     Median value of owner-occupied homes in $1000's\n",
      "\n",
      "    :Missing Attribute Values: None\n",
      "\n",
      "    :Creator: Harrison, D. and Rubinfeld, D.L.\n",
      "\n",
      "This is a copy of UCI ML housing dataset.\n",
      "https://archive.ics.uci.edu/ml/machine-learning-databases/housing/\n",
      "\n",
      "\n",
      "This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.\n",
      "\n",
      "The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic\n",
      "prices and the demand for clean air', J. Environ. Economics & Management,\n",
      "vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics\n",
      "...', Wiley, 1980.   N.B. Various transformations are used in the table on\n",
      "pages 244-261 of the latter.\n",
      "\n",
      "The Boston house-price data has been used in many machine learning papers that address regression\n",
      "problems.   \n",
      "     \n",
      ".. topic:: References\n",
      "\n",
      "   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.\n",
      "   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(boston_data.DESCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sklearn_to_df(sklearn_dataset):\n",
    "    # A helper function to convert the scikit-learn dataset to a pandas DataFrame\n",
    "    # From: https://stackoverflow.com/questions/38105539/how-to-convert-a-scikit-learn-dataset-to-a-pandas-dataset/46379878#46379878\n",
    "    df = pd.DataFrame(sklearn_dataset.data, columns=sklearn_dataset.feature_names)\n",
    "    df['target'] = pd.Series(sklearn_dataset.target)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "    data_url = \"http://lib.stat.cmu.edu/datasets/boston\"\n",
    "    boston = pd.read_csv(data_url, sep=\"\\s+\", skiprows=22, header=None)\n",
    "    #boston = boston.rename(columns={x : i for x, i in enumerate(['CRIM', 'ZN' , 'INDUS' , 'CHAS' , 'NOX' , 'RM' , 'AGE' , 'DIS' , 'RAD' , 'TAX' , 'PTRATIO' , 'B' , 'LSTAT', 'target'])})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.00632</td>\n",
       "      <td>18.00</td>\n",
       "      <td>2.31</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.538</td>\n",
       "      <td>6.575</td>\n",
       "      <td>65.2</td>\n",
       "      <td>4.0900</td>\n",
       "      <td>1.0</td>\n",
       "      <td>296.0</td>\n",
       "      <td>15.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>396.90000</td>\n",
       "      <td>4.98</td>\n",
       "      <td>24.00</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.02731</td>\n",
       "      <td>0.00</td>\n",
       "      <td>7.07</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.469</td>\n",
       "      <td>6.421</td>\n",
       "      <td>78.9</td>\n",
       "      <td>4.9671</td>\n",
       "      <td>2.0</td>\n",
       "      <td>242.0</td>\n",
       "      <td>17.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>396.90000</td>\n",
       "      <td>9.14</td>\n",
       "      <td>21.60</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.02729</td>\n",
       "      <td>0.00</td>\n",
       "      <td>7.07</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.469</td>\n",
       "      <td>7.185</td>\n",
       "      <td>61.1</td>\n",
       "      <td>4.9671</td>\n",
       "      <td>2.0</td>\n",
       "      <td>242.0</td>\n",
       "      <td>17.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          0      1      2    3      4      5     6       7    8      9     10\n",
       "0    0.00632  18.00   2.31  0.0  0.538  6.575  65.2  4.0900  1.0  296.0  15.3\n",
       "1  396.90000   4.98  24.00  NaN    NaN    NaN   NaN     NaN  NaN    NaN   NaN\n",
       "2    0.02731   0.00   7.07  0.0  0.469  6.421  78.9  4.9671  2.0  242.0  17.8\n",
       "3  396.90000   9.14  21.60  NaN    NaN    NaN   NaN     NaN  NaN    NaN   NaN\n",
       "4    0.02729   0.00   7.07  0.0  0.469  7.185  61.1  4.9671  2.0  242.0  17.8"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boston.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "ename": "PatsyError",
     "evalue": "Error evaluating factor: NameError: name 'LSTAT' is not defined\n    target ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT\n                                                                                         ^^^^^",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\compat.py:36\u001b[0m, in \u001b[0;36mcall_and_wrap_exc\u001b[1;34m(msg, origin, f, *args, **kwargs)\u001b[0m\n\u001b[0;32m     35\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m---> 36\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     37\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\eval.py:169\u001b[0m, in \u001b[0;36mEvalEnvironment.eval\u001b[1;34m(self, expr, source_name, inner_namespace)\u001b[0m\n\u001b[0;32m    168\u001b[0m code \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mcompile\u001b[39m(expr, source_name, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124meval\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mflags, \u001b[38;5;28;01mFalse\u001b[39;00m)\n\u001b[1;32m--> 169\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28meval\u001b[39m(code, {}, VarLookupDict([inner_namespace]\n\u001b[0;32m    170\u001b[0m                                     \u001b[38;5;241m+\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_namespaces))\n",
      "File \u001b[1;32m<string>:1\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'LSTAT' is not defined",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mPatsyError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[46], line 2\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mstatsmodels\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mformula\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mapi\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01msmf\u001b[39;00m\n\u001b[1;32m----> 2\u001b[0m est \u001b[38;5;241m=\u001b[39m \u001b[43msmf\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mols\u001b[49m\u001b[43m(\u001b[49m\u001b[43mformula\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mtarget ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[0;32m      3\u001b[0m \u001b[43m              \u001b[49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mboston\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mfit()  \u001b[38;5;66;03m# Does the constant for us\u001b[39;00m\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:200\u001b[0m, in \u001b[0;36mModel.from_formula\u001b[1;34m(cls, formula, data, subset, drop_cols, *args, **kwargs)\u001b[0m\n\u001b[0;32m    197\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m missing \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnone\u001b[39m\u001b[38;5;124m'\u001b[39m:  \u001b[38;5;66;03m# with patsy it's drop or raise. let's raise.\u001b[39;00m\n\u001b[0;32m    198\u001b[0m     missing \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mraise\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m--> 200\u001b[0m tmp \u001b[38;5;241m=\u001b[39m \u001b[43mhandle_formula_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mformula\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdepth\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43meval_env\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    201\u001b[0m \u001b[43m                          \u001b[49m\u001b[43mmissing\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmissing\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    202\u001b[0m ((endog, exog), missing_idx, design_info) \u001b[38;5;241m=\u001b[39m tmp\n\u001b[0;32m    203\u001b[0m max_endog \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mcls\u001b[39m\u001b[38;5;241m.\u001b[39m_formula_max_endog\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\statsmodels\\formula\\formulatools.py:63\u001b[0m, in \u001b[0;36mhandle_formula_data\u001b[1;34m(Y, X, formula, depth, missing)\u001b[0m\n\u001b[0;32m     61\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m     62\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m data_util\u001b[38;5;241m.\u001b[39m_is_using_pandas(Y, \u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m---> 63\u001b[0m         result \u001b[38;5;241m=\u001b[39m \u001b[43mdmatrices\u001b[49m\u001b[43m(\u001b[49m\u001b[43mformula\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mY\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdepth\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreturn_type\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mdataframe\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[0;32m     64\u001b[0m \u001b[43m                           \u001b[49m\u001b[43mNA_action\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mna_action\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     65\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m     66\u001b[0m         result \u001b[38;5;241m=\u001b[39m dmatrices(formula, Y, depth, return_type\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mdataframe\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[0;32m     67\u001b[0m                            NA_action\u001b[38;5;241m=\u001b[39mna_action)\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\highlevel.py:309\u001b[0m, in \u001b[0;36mdmatrices\u001b[1;34m(formula_like, data, eval_env, NA_action, return_type)\u001b[0m\n\u001b[0;32m    299\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Construct two design matrices given a formula_like and data.\u001b[39;00m\n\u001b[0;32m    300\u001b[0m \n\u001b[0;32m    301\u001b[0m \u001b[38;5;124;03mThis function is identical to :func:`dmatrix`, except that it requires\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    306\u001b[0m \u001b[38;5;124;03mSee :func:`dmatrix` for details.\u001b[39;00m\n\u001b[0;32m    307\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m    308\u001b[0m eval_env \u001b[38;5;241m=\u001b[39m EvalEnvironment\u001b[38;5;241m.\u001b[39mcapture(eval_env, reference\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m--> 309\u001b[0m (lhs, rhs) \u001b[38;5;241m=\u001b[39m \u001b[43m_do_highlevel_design\u001b[49m\u001b[43m(\u001b[49m\u001b[43mformula_like\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meval_env\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    310\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43mNA_action\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mreturn_type\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    311\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m lhs\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m    312\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m PatsyError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmodel is missing required outcome variables\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\highlevel.py:164\u001b[0m, in \u001b[0;36m_do_highlevel_design\u001b[1;34m(formula_like, data, eval_env, NA_action, return_type)\u001b[0m\n\u001b[0;32m    162\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mdata_iter_maker\u001b[39m():\n\u001b[0;32m    163\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28miter\u001b[39m([data])\n\u001b[1;32m--> 164\u001b[0m design_infos \u001b[38;5;241m=\u001b[39m \u001b[43m_try_incr_builders\u001b[49m\u001b[43m(\u001b[49m\u001b[43mformula_like\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata_iter_maker\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43meval_env\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    165\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43mNA_action\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    166\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m design_infos \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m    167\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m build_design_matrices(design_infos, data,\n\u001b[0;32m    168\u001b[0m                                  NA_action\u001b[38;5;241m=\u001b[39mNA_action,\n\u001b[0;32m    169\u001b[0m                                  return_type\u001b[38;5;241m=\u001b[39mreturn_type)\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\highlevel.py:66\u001b[0m, in \u001b[0;36m_try_incr_builders\u001b[1;34m(formula_like, data_iter_maker, eval_env, NA_action)\u001b[0m\n\u001b[0;32m     64\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(formula_like, ModelDesc):\n\u001b[0;32m     65\u001b[0m     \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(eval_env, EvalEnvironment)\n\u001b[1;32m---> 66\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mdesign_matrix_builders\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[43mformula_like\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlhs_termlist\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m     67\u001b[0m \u001b[43m                                   \u001b[49m\u001b[43mformula_like\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrhs_termlist\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m     68\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43mdata_iter_maker\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m     69\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43meval_env\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m     70\u001b[0m \u001b[43m                                  \u001b[49m\u001b[43mNA_action\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     71\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m     72\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\build.py:693\u001b[0m, in \u001b[0;36mdesign_matrix_builders\u001b[1;34m(termlists, data_iter_maker, eval_env, NA_action)\u001b[0m\n\u001b[0;32m    689\u001b[0m factor_states \u001b[38;5;241m=\u001b[39m _factors_memorize(all_factors, data_iter_maker, eval_env)\n\u001b[0;32m    690\u001b[0m \u001b[38;5;66;03m# Now all the factors have working eval methods, so we can evaluate them\u001b[39;00m\n\u001b[0;32m    691\u001b[0m \u001b[38;5;66;03m# on some data to find out what type of data they return.\u001b[39;00m\n\u001b[0;32m    692\u001b[0m (num_column_counts,\n\u001b[1;32m--> 693\u001b[0m  cat_levels_contrasts) \u001b[38;5;241m=\u001b[39m \u001b[43m_examine_factor_types\u001b[49m\u001b[43m(\u001b[49m\u001b[43mall_factors\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    694\u001b[0m \u001b[43m                                               \u001b[49m\u001b[43mfactor_states\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    695\u001b[0m \u001b[43m                                               \u001b[49m\u001b[43mdata_iter_maker\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    696\u001b[0m \u001b[43m                                               \u001b[49m\u001b[43mNA_action\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    697\u001b[0m \u001b[38;5;66;03m# Now we need the factor infos, which encapsulate the knowledge of\u001b[39;00m\n\u001b[0;32m    698\u001b[0m \u001b[38;5;66;03m# how to turn any given factor into a chunk of data:\u001b[39;00m\n\u001b[0;32m    699\u001b[0m factor_infos \u001b[38;5;241m=\u001b[39m {}\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\build.py:443\u001b[0m, in \u001b[0;36m_examine_factor_types\u001b[1;34m(factors, factor_states, data_iter_maker, NA_action)\u001b[0m\n\u001b[0;32m    441\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m data \u001b[38;5;129;01min\u001b[39;00m data_iter_maker():\n\u001b[0;32m    442\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m factor \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mlist\u001b[39m(examine_needed):\n\u001b[1;32m--> 443\u001b[0m         value \u001b[38;5;241m=\u001b[39m \u001b[43mfactor\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meval\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfactor_states\u001b[49m\u001b[43m[\u001b[49m\u001b[43mfactor\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    444\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m factor \u001b[38;5;129;01min\u001b[39;00m cat_sniffers \u001b[38;5;129;01mor\u001b[39;00m guess_categorical(value):\n\u001b[0;32m    445\u001b[0m             \u001b[38;5;28;01mif\u001b[39;00m factor \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m cat_sniffers:\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\eval.py:568\u001b[0m, in \u001b[0;36mEvalFactor.eval\u001b[1;34m(self, memorize_state, data)\u001b[0m\n\u001b[0;32m    567\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21meval\u001b[39m(\u001b[38;5;28mself\u001b[39m, memorize_state, data):\n\u001b[1;32m--> 568\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_eval\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmemorize_state\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43meval_code\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    569\u001b[0m \u001b[43m                      \u001b[49m\u001b[43mmemorize_state\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    570\u001b[0m \u001b[43m                      \u001b[49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\eval.py:551\u001b[0m, in \u001b[0;36mEvalFactor._eval\u001b[1;34m(self, code, memorize_state, data)\u001b[0m\n\u001b[0;32m    549\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_eval\u001b[39m(\u001b[38;5;28mself\u001b[39m, code, memorize_state, data):\n\u001b[0;32m    550\u001b[0m     inner_namespace \u001b[38;5;241m=\u001b[39m VarLookupDict([data, memorize_state[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mtransforms\u001b[39m\u001b[38;5;124m\"\u001b[39m]])\n\u001b[1;32m--> 551\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mcall_and_wrap_exc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mError evaluating factor\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[0;32m    552\u001b[0m \u001b[43m                             \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[0;32m    553\u001b[0m \u001b[43m                             \u001b[49m\u001b[43mmemorize_state\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43meval_env\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43meval\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    554\u001b[0m \u001b[43m                             \u001b[49m\u001b[43mcode\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    555\u001b[0m \u001b[43m                             \u001b[49m\u001b[43minner_namespace\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minner_namespace\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mc:\\Python311\\Lib\\site-packages\\patsy\\compat.py:43\u001b[0m, in \u001b[0;36mcall_and_wrap_exc\u001b[1;34m(msg, origin, f, *args, **kwargs)\u001b[0m\n\u001b[0;32m     39\u001b[0m     new_exc \u001b[38;5;241m=\u001b[39m PatsyError(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m     40\u001b[0m                          \u001b[38;5;241m%\u001b[39m (msg, e\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__class__\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, e),\n\u001b[0;32m     41\u001b[0m                          origin)\n\u001b[0;32m     42\u001b[0m     \u001b[38;5;66;03m# Use 'exec' to hide this syntax from the Python 2 parser:\u001b[39;00m\n\u001b[1;32m---> 43\u001b[0m     exec(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mraise new_exc from e\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m     44\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m     45\u001b[0m     \u001b[38;5;66;03m# In python 2, we just let the original exception escape -- better\u001b[39;00m\n\u001b[0;32m     46\u001b[0m     \u001b[38;5;66;03m# than destroying the traceback. But if it's a PatsyError, we can\u001b[39;00m\n\u001b[0;32m     47\u001b[0m     \u001b[38;5;66;03m# at least set the origin properly.\u001b[39;00m\n\u001b[0;32m     48\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(e, PatsyError):\n",
      "File \u001b[1;32m<string>:1\u001b[0m\n",
      "\u001b[1;31mPatsyError\u001b[0m: Error evaluating factor: NameError: name 'LSTAT' is not defined\n    target ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT\n                                                                                         ^^^^^"
     ]
    }
   ],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "est = smf.ols(formula='target ~ CRIM + ZN + INDUS + CHAS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO + B + LSTAT', \n",
    "              data=boston).fit()  # Does the constant for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'est' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[47], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m \u001b[43mest\u001b[49m\u001b[38;5;241m.\u001b[39msummary()\n",
      "\u001b[1;31mNameError\u001b[0m: name 'est' is not defined"
     ]
    }
   ],
   "source": [
    "est.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above table, there is a coef column, which gives the values for $\\beta$ in our model for each independent variable.\n",
    "If the coefficient is negative, there is an inverse relationship between the independent variable and the dependent one.\n",
    "It is important to note that this is not a direct relationship, as retraining the model with just one parameter will likely change this coefficient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "est_simple = smf.ols(formula='target ~ CRIM', \n",
    "              data=boston).fit()  # Does the constant for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Intercept    24.033106\n",
       "CRIM         -0.415190\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est_simple.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the coefficient itself, we are given the standard error, the probability (using the t-statistic) that this value is significant (i.e. if it is less than 0.05), and the lower and upper bounds for the 95% confidence interval - where we can say with 95% confidence that the true value lies within those bounds."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A key reason for this is related to the second warning, indicating there is a strong multicollinearity. We will review this term in the next module and fix the problem it is causing over the next few modules. For now, it indicates that the independent variables are effectively correlated to a high degree, which breaks an assumption with OLS. In short, it means the independent variables are each predicting the same components of the output, and the coefficients are effectively arbitrary. \n",
    "\n",
    "As an example, if we have two variables $a$ and $b$ that are correlated, the coefficient value for $a$ and $b$ in a trained model is effectively shared between them, and whatever value actually appears in the OLS model is just one of many possibilities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the test statistics, good values (for various definitions of \"good\") for these scores allow us to say with a high confidence that the model accurately predicts the data. Bad values indicate that the model should not be used. We will now review a few key values from this table, as a means to validate our model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The $R^2$ statistic\n",
    "\n",
    "The key statistic to review, and the \"one value\" that you are likely to report in your executive summary, is the $R^2$ statistic. It measures how much of the variance in the predicted variable ($Y$) is explained by your model ($X\\beta$), compared to the error of the model ($u$). A high value (near 1) indicates that the model perfectly explains the variable being predicted. A low value (near 0) indicates that the model does not explain the variable at all, which is achieved if the model always predicts the expected value of $Y$. The score can be negative as well, as the model itself can be a net-negative in predictive power (i.e. it model actually predicts incorrectly more than correctly).\n",
    "\n",
    "In the above results, the $R^2$ value is around 0.741, indicating that around 74% of the variance in the predicted variable $Y$ can be explained by the model $X\\beta$. That said, our model has a few problems which we will address soon.\n",
    "\n",
    "To obtain the $R^2$ value, store the regression results object obtained above and extract it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7406077428649428"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.rsquared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "\n",
    "1. Review the documentation at the following link to see what other values can be obtained from a trained estimator:     http://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.RegressionResults.html#statsmodels.regression.linear_model.RegressionResults\n",
    "2. What is the difference between `est.rsquared` and `est.rsquared_adj`? When should you use one over the other?\n",
    "\n",
    "There are quite a few terms on the documentation page we haven't seen yet - many will be reviewed in later modules in this course."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The adjusted R-squared is a modified version of R-squared that accounts for predictors that are not significant in a regression model. In other words, the adjusted R-squared shows whether adding additional predictors improve a regression model or not.\n",
    "\n",
    " - from https://corporatefinanceinstitute.com/resources/knowledge/other/adjusted-r-squared/\n",
    " \n",
    "In general I would prefer the adjusted R-squared unless there is a very simple multiple with a small number of predictors which are uncorrelated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The $F$ statistic\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $F$ statistic is another measure of how significant the fit is. It divides the mean squared error of the model, by the mean squared error of the error term in the model. The probability value under it indicates the probability that we would achieve such a statistic, *if all the coefficients were zero*. In our model, our probability is very low (6.72e-135) indicating there is almost no chance that such an F statistic would be obtained by such a \"zero\" model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "108.07666617432622"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.fvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.722174750114365e-135"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est.f_pvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To put this formally, the F statistic is a test against the null hypothesis:\n",
    "\n",
    "$H_0: \\beta_i = 0 \\forall i$\n",
    "\n",
    "The alternative hypothesis is that *at least* one of the values in $\\beta$ is not 0.\n",
    "\n",
    "The F statistic can be computed using the following terms:\n",
    "\n",
    "$ F = \\frac{ESS}{RSS}$\n",
    "\n",
    "Where $ESS$ is the explained variance of the model and $RSS$ is the unexplained variance. Given the explained variance of the model is due to the component $\\beta X$ and the unexplained component is due to $u$, we can derive the equations as below:\n",
    "\n",
    "$ ESS = \\frac{1}{k-1}\\sum{(\\hat{Y_i} - \\bar{Y})^2}$\n",
    "\n",
    "Where $\\hat{Y_i}$ is the *ith* predicted value and $\\bar{Y}$ is the overall mean of $Y$, and $k$ is the number of variables. In other words, it is the total deviation from the mean that the model explains.\n",
    "\n",
    "For the variance explained by the residuals, we get:\n",
    "\n",
    "$ RSS = \\frac{n}{k}\\sum{u_i^2}$\n",
    "\n",
    "Where $u$ is the error term in our linear regression model and $n$ is the number of samples. There are a few ways to alter these equations to make them easier to compute, all based on performing algebra with the OLS estimator equations defined in earlier modules."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Likelihood Function, Akaike information criterion (AIC) and  Bayesian information criterion (BIC)\n",
    "\n",
    "These three measures are related, and represent the plausibility of the given data given the set of parameters in the model.\n",
    "In all three cases, we use them as relative values. That is, we use these values to compare two different models, and choose the model with the lowest score of these three values (or whichever single statistic you are most concerned with).\n",
    "\n",
    "For instance, if model 1 has a BIC of 3085 and model 2 has a BIC of 4000, we choose model 1.\n",
    "\n",
    "The key function here is the likelihood function, which is used to compute the AIC and BIC. The likelihood function $\\mathcal{L}(\\beta \\mid x)$ is the likelihood that the data could be generated from a model with the given parameters. From an information theory perspective, we aim to maximise the likelihood function. From a computing perspective, it is often easier to both compute the *log likelihood*, and to *minimise the negative log likelihood*. A key component of this is that computers find adding numbers easier than multiplying small numbers, and we can convert from log-space to normal-space using the following pattern:\n",
    "\n",
    "$log(xy) = log(x) + log(y)$\n",
    "\n",
    "When dealing with probabilities, many probability values are very small, and multiplying small numbers near zero is hard for computers. Often, they \"underflow\" and consider a very small number to just be zero, and then any product from that point on is zero. Instead, we compute the log of all numbers and add them together - no underflow!\n",
    "\n",
    "Once the likelihood function (or negative log likelihood) has been computed, the maximum value it can take (when optimised) is $\\hat{L}$. From here, the AIC is defined as:\n",
    "\n",
    "$ AIC = sk - s\\ln(\\hat{L})$\n",
    "\n",
    "The BIC is defined similarly:\n",
    "\n",
    "$ BIC = \\ln(n)k - 2\\ln({\\hat L})$\n",
    "\n",
    "Typically the BIC is preferred, as it is more stable in most circumstances. However, for the BIC to be valid, the number of samples must be much more than the number of parameters."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
